c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      SUBROUTINE USRMAT(RZE,DRZE,NONP,NDIM,NSTS,NEL,IPT, DSTRAN,DFGR0, 
     +    DFGR1,ROT1,STRESS,EPX,TOLDG,TNEWG,TREF,DSAVE,TIM,DTIME,USRVE,  
     +    NUSRVE,USRVN,NUSRVN,VELGRD,IMAT,KSTP,ITR,AXI,   
     +    EFSTS,YPS,FIP,IND)  
C======================================================================
C
C     User constitutive law subroutine - UMAT
C
C======================================================================
C     
C     Material Properties
C     can be accessed using user-defined variable 
C
C---------------
C     INPUT
C---------------
C     RZE            : Nodal Point Coordinates   
C     DRZE          : Nodal Point Displacements
C     NONP          : Number of Nodes in The Element  
C     NDIM           : Space Dimension      
C     NSTS           :  Size of the Stress or Strain Component Array         
C     NEL             : Element Number
C     IPT              : Integration Point Number     
C     DSTRAN      : Total strain increment from n to n+1 (possitioned at t=n)
C     DFGR0         : Deformation gradient F(t=n) w.r.t. t=0  
C     DFGR1         : Deformation gradient F(t=n+1) w.r.t. t=0
C     ROT1          : Rotation tensor from n to n+1/2
C     STRESS      : Stress at the start of increment (possitioned at t=n)  
C     EPX            : Effective plastic strain at the start of increment 
C     TOLDG       : Temperature at the previous increment
C     TNEWG      : Temperature at the current increment
C     TREF         : Reference temperature 
C     TIM(1)       : 0.0
C     TIM(2)       : Current total time (i.e. at t=n)
C     DTIME       : Time increment
C     USRVE       : User-Defined Elemental Variables at the start of increment  
C     NUSRVE     : Number of User-Defined Elemental Variables
C     USRVN       :User-Defined Nodal Variables
C     VELGRD     : Velocity Gradients dv/dx 
C     NUSRVN     : Number of User-Defined Nodal Variables
C     IMAT          :Material Group Number (Used for Multiple User Materials,
C                           Will be Replaced by User Routine Number in The Future)
C     KSTP         : Step Number
C     ITR           : Iteration Number
C     AXI           : Material Axes     
C     IND          : Index which specifies the object type (Only 3D)
C----------------
C     OUTPUT
C----------------
C     STRESS       :  stress at the end increment
C     DSAVE        : Tangent modulus at the end of increment
C     USRVE        : User-defined Elemental Variables at the end of increment  
C     EFSTS         : FLOW STRESS (Only 3D)   
C     YPS:    DERIVATIVE OF FLOW STRESS W.R.T. EFFECTIVE PLASTIC STRAIN (Only 3D)
C     FIP:  DERIVATIVE OF FLOW STRESS W.R.T. EFFECTIVE PLASTIC STRAIN RATE (Only 3D) 
C=====================================================================

      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
      DIMENSION STRESS(NSTS),DSAVE(NSTS,NSTS),DSTRAN(NSTS),TIM(2),   
     +                 DFGR0(NDIM,NDIM),DFGR1(NDIM,NDIM),   
     +                 ROT1(NDIM,NDIM),USRVE(NUSRVE),RZE(NDIM,NONP),  
     +                 DRZE(NDIM,NONP),USRVN(NUSRVN,NONP),   
     +                 VELGRD(NDIM,NDIM),AXI(NDIM,NDIM)   

C    TO BE WRITTEN BY USER
c
      return
      end

C======================================================================
      SUBROUTINE USRMAT_E(RZE,DRZE,NONP,NDIM,NSTS,NEL,IPT,DSTRAN,DFGR0,  
     +    DFGR1,ROT1,STRESS,EPX,TOLDG,TNEWG,TREF,DSAVE,TIM,DTIME,USRVE,  
     +    NUSRVE,USRVN,NUSRVN,VELGRD,IMAT,KSTP,ITR,AXI,   
     +    EFSTS,YPS,FIP,IND) 
C======================================================================
C
C     2D user subroutine example for linear elastic material 
C
C======================================================================
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)

      DIMENSION STRESS(4),DSAVE(4,4),DSTRAN(4),TIM(2),   
     +                 DFGR0(2,2),DFGR1(2,2),ROT1(2,2),USRVE(NUSRVE),   
     +                 RZE(2,4),DRZE(2,4),USRVN(NUSRVN,4),  
     +                 VELGRD(2,2),AXI(2,2)        
      
      YM=1000
      PV=0.3
      THIRD=-1./3.
C
C     Elastic tangent modulus
C     
      S1=(1.+PV)*(1.0-2.0*PV)
      S2=(1.-PV)/S1
      S3=YM/(1.+PV)
      S4=S3/2.
      S5=S3/(1.-2.*PV)
      S6=S5*PV
      S7=S5-S6
C
C     Store elastic tangent modulus in DSAVE
C
      DSAVE(1,1)=S7
      DSAVE(1,2)=S6
      DSAVE(1,3)=S6
      DSAVE(2,1)=S6
      DSAVE(2,2)=S7
      DSAVE(2,3)=S6
      DSAVE(4,4)=S4
      DSAVE(3,1)=S6
      DSAVE(3,2)=S6
      DSAVE(3,3)=S7
      DSAVE(4,1)=0.
      DSAVE(1,4)=0.
      DSAVE(4,2)=0.
      DSAVE(2,4)=0.
      DSAVE(4,3)=0.
      DSAVE(3,4)=0.
C     
C     Total stain increment 
C     
      H1=DSTRAN(1)
      H2=DSTRAN(2)
      H3=DSTRAN(3)
      H4=DSTRAN(4)
C
C     Calulate elastic stress
C     
      T1=STRESS(1)+S7*H1+S6*H2+S6*H3
      T2=STRESS(2)+S6*H1+S7*H2+S6*H3
      T3=STRESS(3)+S6*H1+S6*H2+S7*H3
      T4=STRESS(4)+S4*H4
C
C     new stress at the end of increment
C     
      STRESS(1)=T1
      STRESS(2)=T2
      STRESS(3)=T3
      STRESS(4)=T4
C
C     No effective plastic strain 
C
      EPX=0.0
C
      RETURN
      END
C
C************************************************************
